{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f91d538a-90e1-49f2-8b3f-9777c17b7770",
   "metadata": {},
   "source": [
    "# Classic VRPs\n",
    "\n",
    "This notebook shows how to use PyVRP to solve two classic variants of the VRP: the capacitated vehicle routing problem (CVRP), and the vehicle routing problem with time windows (VRPTW).\n",
    "It builds on the tutorial by solving much larger instances, and going into more detail about the various plotting tools and diagnostics available in PyVRP.\n",
    "\n",
    "A CVRP instance is defined on a complete graph $G=(V,A)$, where $V$ is the vertex set and $A$ is the arc set. \n",
    "The vertex set $V$ is partitioned into $V=\\{0\\} \\cup V_c$, where $0$ represents the depot and $V_c=\\{1, \\dots, n\\}$ denotes the set of $n$ customers.\n",
    "Each arc $(i, j) \\in A$ has a weight $d_{ij} \\ge 0$ that represents the travel distance from $i \\in V$ to $j \\in V$.\n",
    "Each customer $i \\in V_c$ has a demand $q_{i} \\ge 0$.\n",
    "The objective is to find a feasible solution that minimises the total distance.\n",
    "\n",
    "A VRPTW instance additionally incorporates time aspects into the problem.\n",
    "For the sake of exposition we assume the travel duration $t_{ij} \\ge 0$ is equal to the travel distance $d_{ij}$ in this notebook.\n",
    "Each customer $i \\in V_c$ has a service time $s_{i} \\ge 0$ and a (hard) time window $\\left[e_i, l_i\\right]$ that denotes the earliest and latest time that service can start.\n",
    "A vehicle is allowed to arrive at a customer location before the beginning of the time window, but it must wait for the window to open to start the delivery. \n",
    "Each vehicle must return to the depot before the end of the depot time window $H$.\n",
    "The objective is to find a feasible solution that minimises the total distance.\n",
    "\n",
    "Let's first import what we will use in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84952b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from tabulate import tabulate\n",
    "from vrplib import read_solution\n",
    "\n",
    "from pyvrp import Model, read\n",
    "from pyvrp.plotting import (\n",
    "    plot_coordinates,\n",
    "    plot_instance,\n",
    "    plot_result,\n",
    "    plot_route_schedule,\n",
    ")\n",
    "from pyvrp.stop import MaxIterations, MaxRuntime"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "de3370c7",
   "metadata": {},
   "source": [
    "## The capacitated VRP"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d48d3c9d",
   "metadata": {},
   "source": [
    "### Reading the instance\n",
    "\n",
    "We will solve the `X-n439-k37` instance, which is part of the [X instance set](http://vrp.atd-lab.inf.puc-rio.br/index.php/en/new-instances) that is widely used to benchmark CVRP algorithms.\n",
    "The function `pyvrp.read` reads the instance file and converts it to a `ProblemData` instance. \n",
    "We pass the argument `round_func=\"round\"` to compute the Euclidean distances rounded to the nearest integral, which is the convention for the X benchmark set.\n",
    "We also load the best known solution to evaluate our solver later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc6a35a-3cee-4371-b62a-c80f044cf07f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "INSTANCE = read(\"data/X-n439-k37.vrp\", round_func=\"round\")\n",
    "BKS = read_solution(\"data/X-n439-k37.sol\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d0a3bb5e-01ea-4738-b327-c4aec609e37a",
   "metadata": {},
   "source": [
    "Let's plot the instance and see what we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67f7ef22-7925-4e32-8f2a-59d92af0a4a3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(figsize=(8, 8))\n",
    "plot_coordinates(INSTANCE, ax=ax)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "71ebf6cc-a34f-46ef-9606-85d86b2c0ef3",
   "metadata": {},
   "source": [
    "### Solving the instance\n",
    "\n",
    "We will again use the `Model` interface to solve the instance.\n",
    "The `Model` interface supports a convenient `from_data` method that can be used to instantiate a model from a known `ProblemData` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b7883a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model.from_data(INSTANCE)\n",
    "result = model.solve(stop=MaxIterations(2000), seed=42, display=False)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5179830f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gap = 100 * (result.cost() - BKS[\"cost\"]) / BKS[\"cost\"]\n",
    "print(f\"Found a solution with cost: {result.cost()}.\")\n",
    "print(f\"This is {gap:.1f}% worse than the best known\", end=\" \")\n",
    "print(f\"solution, which is {BKS['cost']}.\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d23a2374-15f5-4e8b-b021-a6c8eae1c619",
   "metadata": {},
   "source": [
    "We've managed to find a very good solution quickly!\n",
    "\n",
    "The `Result` object also contains useful statistics about the optimisation.\n",
    "We can now plot these statistics as well as the final solution use `plot_result`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5bdaff4-adbb-4e58-a401-b56b8563ff19",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 9))\n",
    "plot_result(result, INSTANCE, fig)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "250b1595",
   "metadata": {},
   "source": [
    "PyVRP internally uses a genetic algorithm consisting of a population of feasible and infeasible solutions.\n",
    "These solutions are iteratively combined into new offspring solutions, that should result in increasingly better solutions. \n",
    "Of course, the solutions should not all be too similar: then there is little to gain from combining different solutions.\n",
    "The top-left *Diversity* plot tracks the average diversity of solutions in each of the feasible and infeasible solution populations.\n",
    "The *Objectives* plot gives an overview of the best and average solution quality in the current population.\n",
    "The bottom-left figure shows iteration runtimes in seconds.\n",
    "Finally, the *Solution* plot shows the best observed solution."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4846b4b8",
   "metadata": {},
   "source": [
    "## The VRP with time windows\n",
    "\n",
    "### Reading the instance\n",
    "\n",
    "We start with a basic example that loads an instance and solves it using the standard configuration used by the `Model` interface.\n",
    "For the basic example we use one of the well-known Solomon instances.\n",
    "\n",
    "We again use the function `pyvrp.read`. We pass the argument `round_func=\"dimacs\"` following the [DIMACS VRP challenge](http://dimacs.rutgers.edu/programs/challenge/vrp/) convention, this computes distances and durations truncated to one decimal place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c847d21",
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE = read(\"data/RC208.vrp\", round_func=\"dimacs\")\n",
    "BKS = read_solution(\"data/RC208.sol\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "fe1add17",
   "metadata": {},
   "source": [
    "Let's plot the instance and see what we have.\n",
    "The function `plot_instance` will plot time windows, delivery demands and coordinates, which should give us a good impression of what the instance looks like.\n",
    "These plots can also be produced separately by calling the appropriate `plot_*` function: see [the API documentation](https://pyvrp.readthedocs.io/en/latest/api/plotting.html) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bb7a41a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 6))\n",
    "plot_instance(INSTANCE, fig)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cabf0ded",
   "metadata": {},
   "source": [
    "### Solving the instance\n",
    "\n",
    "We will again use the `Model` interface to solve the instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66db0c08",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model.from_data(INSTANCE)\n",
    "result = model.solve(stop=MaxIterations(1000), seed=42, display=False)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bcc878b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cost = result.cost() / 10\n",
    "gap = 100 * (cost - BKS[\"cost\"]) / BKS[\"cost\"]\n",
    "print(f\"Found a solution with cost: {cost}.\")\n",
    "print(f\"This is {gap:.1f}% worse than the optimal solution,\", end=\" \")\n",
    "print(f\"which is {BKS['cost']}.\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "428ffeb1",
   "metadata": {},
   "source": [
    "We've managed to find a (near) optimal solution in a few seconds!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1503f076",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 9))\n",
    "plot_result(result, INSTANCE, fig)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "731869e3",
   "metadata": {},
   "source": [
    "We can also inspect some statistics of the different routes, such as route distance, various durations, the number of stops and total delivery amount."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dd99e17",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = result.best\n",
    "routes = solution.routes()\n",
    "\n",
    "data = [\n",
    "    {\n",
    "        \"num_stops\": len(route),\n",
    "        \"distance\": route.distance(),\n",
    "        \"service_duration\": route.service_duration(),\n",
    "        \"wait_duration\": route.wait_duration(),\n",
    "        \"time_warp\": route.time_warp(),\n",
    "        \"delivery\": route.delivery(),\n",
    "    }\n",
    "    for route in routes\n",
    "]\n",
    "\n",
    "header = list(data[0].keys())\n",
    "rows = [datum.values() for datum in data]\n",
    "tabulate(rows, header, tablefmt=\"html\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9a3fb318",
   "metadata": {},
   "source": [
    "We can inspect the routes in more detail using the `plot_route_schedule` function.\n",
    "This will plot distance on the x-axis, and time on the y-axis, separating actual travel/driving time from waiting and service time.\n",
    "The clients visited are plotted as grey vertical bars indicating their time windows.\n",
    "We can see a jump to the start of the time window in the main (earliest) time line when a vehicle arrives early at a customer and has to wait.\n",
    "In some cases, there is slack in the route indicated by a semi-transparent region on top of the earliest time line.\n",
    "The grey background indicates the remaining load of the truck during the route, where the (right) y-axis ends at the vehicle capacity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4afa476d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axarr = plt.subplots(2, 2, figsize=(15, 9))\n",
    "for idx, (ax, route) in enumerate(zip(axarr.reshape(-1), routes)):\n",
    "    plot_route_schedule(\n",
    "        INSTANCE,\n",
    "        route,\n",
    "        title=f\"Route {idx}\",\n",
    "        ax=ax,\n",
    "        legend=idx == 0,\n",
    "    )\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8182b4f6",
   "metadata": {},
   "source": [
    "## Solving a larger VRPTW instance\n",
    "\n",
    "To show that PyVRP can also handle much larger instances, we will solve one of the largest Gehring and Homberger VRPTW benchmark instances.\n",
    "The selected instance - `RC2_10_5` - has 1000 clients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db87047f",
   "metadata": {},
   "outputs": [],
   "source": [
    "INSTANCE = read(\"data/RC2_10_5.vrp\", round_func=\"dimacs\")\n",
    "BKS = read_solution(\"data/RC2_10_5.sol\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4522453",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 9))\n",
    "plot_instance(INSTANCE, fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89134da9",
   "metadata": {},
   "source": [
    "Here, we will use a runtime-based stopping criterion: we give the solver 30 seconds to compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa29e516",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model.from_data(INSTANCE)\n",
    "result = model.solve(stop=MaxRuntime(30), seed=42, display=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ba2be7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cost = result.cost() / 10\n",
    "gap = 100 * (cost - BKS[\"cost\"]) / BKS[\"cost\"]\n",
    "print(f\"Found a solution with cost: {cost}.\")\n",
    "print(f\"This is {gap:.1f}% worse than the best-known solution,\", end=\" \")\n",
    "print(f\"which is {BKS['cost']}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "455ccf31",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_result(result, INSTANCE)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "59a499d0-695e-45ab-87fc-2382138bc744",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In this notebook, we used PyVRP's `Model` interface to solve a CVRP instance with 438 clients to near-optimality, as well as several VRPTW instances, including a large 1000 client instance.\n",
    "Moreover, we demonstrated how to use the plotting tools to visualise the instance and statistics collected during the search procedure."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
